arXiv: 1504.04865vl [astro-ph.IM] 19 Apr 2015 


A new method for deriving the stellar birth function 
of resolved stellar populations^ 


M. Gennaro^, K. Tchernyshyov^, T. M. Brown^, K. D. Gordon^’^ 

gennaroSstsci.edu 

^Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218 
^Department of Physics and Astronomy, The Johns Hopkins University, 

3400 N. Charles Street, Baltimore, MD 21218 
^ Sterrenkundig Observatorium, Universiteit Gent, Gent, Belgium 


ABSTRACT 

We present a new method for deriving the stellar birth function (SBF) of resolved stellar 
populations. The SBF (stars born per unit mass, time, and metallicity) is the combination of the 
initial mass function (IMF), the star-formation history (SFH), and the metallicity distribution 
function (MDF). The framework of our analysis is that of Poisson Point Processes (PPPs), a class 
of statistical models suitable when dealing with points (stars) in a multidimensional space (the 
measurement space of multiple photometric bands). The theory of PPPs easily accommodates 
the modeling of measurement errors as well as that of incompleteness. Compared to most of the 
tools used to study resolved stellar populations, our method avoids binning stars in the color- 
magnitude diagram and uses the entirety of the information (i.e., the whole likelihood function) 
for each data point; the proper combination of the individual likelihoods allows the computation 
of the posterior probability for the global population parameters. This includes unknowns such 
as the IMF slope and combination of SFH and MDF, which are rarely solved for simultaneously 
in the literature, however entangled and correlated they might be. Our method also allows 
proper inclusion of nuisance parameters, such as distance and extinction distributions. The aim 
of this paper, is to assess the validity of this new approach under a range of assumptions, using 
only simulated data. Forthcoming work will show applications to real data. Although it has a 
broad scope of possible applications, we have developed this method to study multi-band HST 
observations of the Milky Way Bulge. Therefore we will focus on simulations with characteristics 
similar to those of the Galactic Bulge . 

Subject headings: methods: statistical, Galaxy: bulge 


1. Introduction 

The study of resolved stellar populations is go¬ 
ing through a remarkable growth period, with 
space observatories like the Hubble Space Tele¬ 
scope (HST) providing high-resolution probes 
through nearby galaxies, and all-sky surveys 
like the Sloan Digital Sky Survey (SDSS) and 
Panoramic Survey Telescope and Rapid Response 


^Based on observations made with the NASA/ESA Hub¬ 
ble Space Telescope, obtained at STScI, which is operated 
by AURA, Inc., under NASA contract NAS 5-26555 


System (Pan-STARRS) revealing ever more sub¬ 
structure in the Local Group. Upcoming missions 
like the James Webb Space Telescope (JWST) and 
the Large Synoptic Survey Telescope (LSST) will 
further increase this wealth of data. A big chal¬ 
lenge for the current and next generations of as¬ 
tronomers is that of developing appropriate tools 
to make the best out of these data, with sound 
statistical methods to help understand the un¬ 
derlying errors and properly interpret the results. 
Ultimately, in the context of stellar populations, 
this means correctly interpreting the features ob- 
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served in the color-magnitude diagram (CMD), or 

nandez et al. 1999 

jjprgensen & Lindegren 

2005 

in equivalent diagrams. 

Naylor & Jeffries 

2006 Da Rio et al.|2010[ 

Walm-| 

The recent literature offers several examples of 

swell et al. 2013 

). The last paper also gives an 


such methods that can be grouped into two broad 
categories: methods based on binning the data us¬ 
ing a grid defined on the CMD, and methods that 
do not bin the data, but instead use each individ¬ 
ual measurements separately. The latter, at least 
in some respects, can be seen as the limit of the 
former in the case of very small grid cells where 
only 0 or 1 stars are observed in each cell. 

Most of the bin-based methods, with their own 
differences, follow a similar approach. They first 
create a set of basis functions, i.e., simulated 
CMDs of simple stellar populations; these are usu¬ 
ally realized with Monte Carlo techniques, ac¬ 
counting for the photometric errors and selection 
effects. They then linearly combine the basis func¬ 
tions to produce a synthetic CMD for the whole 
population. The fit is then performed via mini¬ 
mization of statistics comparing the predicted and 
observed numbers of stars in each grid-cell, in or¬ 
der to find the appropriate weight for each ba¬ 
sis function. These weights correspond to the in¬ 
tensity of the star formation episode associated 
with each simple stellar population. Examples of 


bin-based methods are 

those developed by|Harris 

& Zaritsky (2001 

), Vergely et al. (2002), Dolphin 

(2002), Ng et al. 

2002 

, Cignoni et al. (2006), and 

Aparicio & Hidalgo (2 

009). These methods have 


been successfully applied to a wide range of ob¬ 
servations. However, they become harder to ap¬ 
ply when very few stars are observed, making a 
large number of CMD cells empty, or equivalently, 
forcing the use of very large cells. This limitation 
is mostly evident when more than 2 photometric 
bands data are available, making the number of 
useful cells even smaller, a problem that can be 
seen as one form of the curse of dimensionality. 

Unbinned methods try and use the full infor¬ 
mation available for each datum, taking into ac¬ 
count the noise associated with individual mea¬ 
surements. Generally speaking, the methods in 
this class are based on computing the probability 
of each observed datum given the available stel¬ 
lar evolutionary models. The individual probabil¬ 
ities are then appropriately combined to derive the 
population parameters. Again, several such meth¬ 
ods are described in the literature, each with its 


insightful and detailed review of the differences be¬ 
tween (and within) binned and unbinned methods. 

To build our method we start by recognizing 
that, on a population scale, the star formation 
process can be described as a stochastic process 
in which individual stars are drawn independently 
from a parent distribution; this kind of process can 
be modeled as a type of random process known as 
a Poisson Point Process (PPP). If we define our 
PPP on the space of intrinsic physical parameters, 
i.e. if we consider stellar age, mass, and metallic- 
ity as the stochastic variables, the parent distri¬ 
bution is simply the stellar birth function (SBF), 
which represents the probability of there being 
a star with a certain age, mass, and metallicity. 
The SBF can be regarded as a combination of the 
star formation history (SFH), initial mass func¬ 
tion (IMF), and metallicity distribution function 
(MDF). 

The key to our method is that in the PPP for¬ 
malism, one can map the probability distribution 
on the intrinsic parameter space to an equivalent 
probability distribution on a space of arbitrary ob¬ 
servable quantities, such as (noisy) photometric 
observations or individual stellar spectra. Addi¬ 
tionally, this mapping can include the fact that 
some stars, while born and hence relevant to the 
SBF, cannot be observed; this part of the map¬ 
ping is related to data incompleteness and is tech¬ 
nically referred to as thinning. For example, these 
stars may have evolved off the main sequence or 
be fainter than the detection limit of the specific 
observations. The form of this mapping is such 
that for each observed star, we can easily com¬ 
pute the likelihood of any combination of intrinsic 
parameters, including nuisance parameters such as 
distance or reddening. These likelihoods are then 
combined to compute the posterior probability of 
any given SBF in a way that fully uses each indi¬ 
vidual measurement’s information. 

The statistical model we derive is similar to 


the one by Weisz et al. (2013), where the au¬ 


own peculiarities (e.g. Tolstoy & Saha 1996 Her- 


thors use a hierarchical Bayesian approach to ob¬ 
tain the posterior probability of the slope of the 
high-mass IMF in the context of stellar clusters 
analysis. Other than using a different path for 
reaching a similar model description, there are fur- 
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ther differences between our approach and that by 


Weisz et al. (2013). For example, we are not be¬ 


ing limited to studying a single parameter, since 
we include the whole SFH and MDF in addition 
to the IMF slope. Moreover, we simulated data to 
test our method, taking into account both noise 
and incompleteness and we have developed and 
described a full numerical approach for the actual 
calculation of both the likelihood and complete¬ 
ness functions. This is in contrast to the approach 
of Weisz et al. (2013) who do not derive the like¬ 
lihoods from photometric measurements, but in¬ 
stead use an analytic approximation to describe 
the likelihoods of individual objects as well as the 
incompleteness function. Within the astronomical 
community PPPs have been used to study prob¬ 
lems in other areas; without the presumption of 
being exhaustive, but only to show the broad scope 


Tremaine ( 

2002) 

Youdin 

(2011) 

Foreman-Mackey 

et al. 

(201 

4) in the area of exoplanet search and 


population census, Lombardi et al. (2013) who use 
PPPs to study the local Schmidt law in molecular 
clouds, and Hugeback et al. (2007) who model the 
Quasar Luminosity function in magnitude-redshift 
space as a PPP. 


Our implementation of PPPs for SBF deter¬ 
mination requires a discretization of the intrin¬ 
sic physical parameter space on a grid. This is 
common to all the existing methods since all the 
mappings between physical parameters and obser¬ 
vations are based on stellar models, which only 
exist for finite grids of parameters. Our method 
formally accounts for such approximation in its 
definition. It is important to note that if the 
adopted grid of models is fine enough to resolve 
each star’s likelihood function, the discretized ap¬ 
proach is substantially equivalent to a completely 
continuous one. 


The PPP formalism has several important ad¬ 
vantages: (i) it is an exact and faithful mathe¬ 
matical analogue to the generally accepted idea 
of stellar population formation; (ii) it allows us 
to exploit synergies between subfields of astron¬ 
omy and between astronomy and applied mathe¬ 
matics and statistics; (iii) it is conveniently mod¬ 
ular; and (iv) it is flexible and extendible. Ex¬ 
actness (i) makes the formalism a useful and as- 
trophysically motivated starting point for devel¬ 
oping practical techniques. In particular, all of 


the existing CMD fitting techniques mentioned 
above can be derived from the fully general PPP 
formalism by taking various combinations of sim¬ 
plifying assumptions, approximations, and limits. 
Intra- and inter-disciplinary synergy (ii) simplifies 
the process of devising new computational tech¬ 
niques and verifying old ones. For example, we 
use lessons learned from medical imaging, specif¬ 
ically positron emission tomography to find the 
best-fit SBF. Other investigators can now apply 
the vast statistical literature on different types 
of optimal approximate methods to vet existing 
methods (e.g. to show that they are theoretically 
unbiased and have good variance properties) and 
develop new ones. As we will show in the paper, 
the modularity of the SBF posterior probability 
is computationally convenient (iii), since it allows 
us to separately precompute several of the nec¬ 
essary quantities (e.g., the individual likelihoods) 
and makes the other computations involved paral- 
lelizable and therefore fast. This modularity also 
makes the formalism and method flexible and ex¬ 
tendable (iv). For example, we use a specific set 
of stellar evolutionary libraries and photometric 
bands, but the formalism and method both ap¬ 
ply to any set of libraries and photometric bands. 
Moreover, the method can also be expanded to 
include other kinds of observations, such as spec¬ 
tra of individual stars, purely by modifying the 
likelihood and thinning functions. We will show 
examples in which spectroscopic constraints are 
incorporated to provide strong information on the 
populations’ MDF. 


In future work, we intend to exploit our 
method’s flexibility to analyze the data from the 


Galactic Bulge Treasury Program (Brown et al. 


2009, 2010). This is a deep HST dataset which in¬ 


cludes five photometric bands, is supplemented by 
available ground-based stellar spectroscopy, and 
which targets four fields in the Galactic Bulge, 
where there exists significant star-to-star distance 
and reddening variation. 


The paper is structured as follows: we describe 
the basics of Poisson Point Processes in Section [2] 
We then describe the adopted library of stellar 
models in Section |3l Section |4] deals with the 
treatment of measurement errors and incomplete¬ 
ness and how they affect the specifications of the 
individual likelihoods. Section [S] describes the ex¬ 
plicit solution of the population properties. We 
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outline the test catalogs simulation process in Sec¬ 
tion and apply our method on such catalogs and 
show the results in Section [T] We summarize our 
findings in Section 


2. Poisson Point Processes 


A Poisson Point Process (PPP) is a statisti¬ 
cal model that describes the counting of points 
in a multi-dimensional space. A full description 
of PPPs and some of their applications is given in 


Streit (2010). Given the unfamiliar nature of these 


models to the astronomical community, in the fol¬ 
lowing we summarize the basics of PPPs following 
the development in Streit (20101. In particular, 
we focus on the aspects of PPPs that are relevant 
to the analysis of stellar populations. 


A realization of a PPP consists of a certain 
number of points observed in a state space S. For 
our purposes, the state space is a 3D space in 
which the coordinates are stellar mass, age, and 
metallicitjQ We will show in the following how the 
PPP of interest in the study of the SBF, which is 
defined in the stellar parameters space, can be re¬ 
lated to a different PPP, i.e. the noisy incomplete 
set of photometric measurement we have access to. 


One can easily imagine extending S to include 
other dimensions that can be of interest in the 
study of stellar populations, such as stellar dis¬ 
tance, reddening, and multiplicity properties. The 
formalism would be equivalent and, for the sake 
of simplicity, we will drop these additional dimen¬ 
sions. While the state space S can in principle be 
infinite, we are interested in realizations of PPPs 
on a bounded subset TZ of S. For the purpose of 
stellar population analysis, the subset TZ will be 
defined by the range of stellar parameters in the 
specific set of adopted models (see Sect.[^. 

In describing a PPP, both the number and the 
distribution of points over the state space are ran¬ 
dom variables. A realization of a PPP in TZ is 


^Even if we denote the metallicity with the lowercase sym¬ 
bol 2 , we always make use of the spectroscopic definition 
of metallicity: [M/H] = log ^ ^ — log ^ ^ ^ . 

Because models can include variations in o-element abun¬ 
dances, [Fe/H] and [M/H] are not necessarily the same. 
However, because a given set of models is computed us¬ 
ing a well-defined [a/Fe] vs. [Fe/H] relation, [Fe/H] will be 
sometimes used or mentioned instead of [M/H]. 


denoted by 


? = • • • ,Sn}) (1) 

where the total number of points n is explicitly 
indicated and the ordering of the Si points is irrel¬ 
evant; the Si can, in principle, include duplicates. 
The most important quantity that characterizes a 
PPP is its intensity, A(s), which describes how the 
Si points are distributed in the state space. The in¬ 
tensity must be non-negative everywhere and has 
to integrate to a finite value over the state space: 


Vs S S, A(s) >0; 0 < / A(s)ds < oo (2) 

Jn 


Two conditions must be met for a model to be 
a PPP: 1) the total number of points in the subset 
TZ, N-ji, has a Poisson distribution with parameter 
A(s)ds and 2) if TZi and TZ 2 are disjoint, then 
and /Vt^^ independent. 

Having set the stage, we show how a realization 
^ can be generated. First, the number of points n 
is a draw of the Poisson variable, N, distributed 
according to: 

PN{n) = ; p = [ A(s)ds . (3) 

n'- Jn 

The integral of the intensity therefore gives the 
expected number of points. The location of these 
points in TZ is given by n independent draws of the 
random variable S with probability distribution 
function (pdf) given by 


Ps{s) 


Ks) 

In 


( 4 ) 


We introduce a random variable S = {N,X), 
where N is the number of points and X = 
{si,--- ,s„} is the points set. The probability 
of a generic event evaluated at S = is given by 

Ph(0 =PN{n)px\N{{si,--- ,Sn}\n) (5) 

where the first factor is given by Eq. ([^, the sec¬ 
ond is 

n 

Px\Ni {si, • • • ,Sn}\n) = n! ]^ps(si): (6) 

and ps is the pdf for a single point, as given by 
Eq. The n! factor is due to the fact that there 
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are n! possible combinations of the s^’s that corre¬ 
spond to the unordered set X. Combining every¬ 
thing we obtain that 

n 

Z=1 

In the context of stellar populations, the inten¬ 
sity function corresponds to the stellar birth func¬ 
tion, i.e., represents how many stars have been 
formed per unit time, mass, and metallicity, within 
the patch of the sky under study. This function is 
a combination of the initial mass function (IMF), 
the star-formation history (SFH), and the metal¬ 
licity distribution function (MDF). While in prin¬ 
ciple we might expect the IMF to be a function 
of both age and metallicity, in the following we 
will make the simplifying assumption that it is in¬ 
stead independent with respect to both. Therefore 
A(s) = A(m, a, z) oc I{m) $(a, z). 

Equation Q gives the probability of n stars 
being formed with a given set of properties 
{(m, a, z)i, • • • ,{m,a,z)n}, i-e., it is the proba¬ 
bility of the set of physical parameters {m,a,z), 
given A; pe{^) = ps{^\X)- Using Bayes’ theo¬ 
rem, we can write the probability of A given the 
available {{m,a, z)}: 

Pa{X\0 (xpe{^\X)pa{X) ( 8 ) 

where the normalization factor is omitted. If we 
had the parameters (to, a, z) for all the stars in our 
sample, our computational problem would be ex¬ 
ploring the pdf of A defined in Eq. ([^. However, 
the problem to solve in the case of stellar popu¬ 
lations is that of determining A given the number 
of observed stars and a set of flux or magnitude 
measurements for each star. The physical param¬ 
eters (to, a, z) are not directly observable or acces¬ 
sible, so stellar models must be used to interpret 
the measured quantities in terms of (to, a, z). The 
uncertainties in these physical parameters will de¬ 
pend upon the measurement and modeling errors. 
Furthermore, of the stars that are born according 
to the true, underlying A, only some will be ob¬ 
servable. The incompleteness of the data is due to 
both stellar evolution (massive stars ending their 
lives as stellar remnants) and observational limits 
(the ability to detect an object and of measure its 
flux given the noise and crowding properties of the 
specific observations). 


The theory of PPPs, as outlined in Streit 


( 20101 , easily accommodates measurement pro¬ 
cesses and incompletenes^ We summarize the 
treatment of both below. 


2.1. Measurement process and errors 

The unknown intensity A assumes values in 
the space S of stellar mass, age and metallicity. 
The actual measurements are made in another 
space, which we indicate with T and which can be 
thought of as a A:-dimensional magnitude space, 
where k is the number of available bands. 0 

For a set of stellar parameters s = (to, a, z), we 
can determine the probability p{ t\ s) of observing 
a set of magnitudes t = (Mi, • • • , Mk) using stel¬ 
lar models and observational uncertainties. Given 
a realization ^ = (n, {(to, a, z)i, • • • , {m,a, z)n}) 
of a PPP with intensity A, it can be shown that 
77 = (n, {(Ml, • • • , Mfe)i, • • • , (Ml, • • • , Mk)n}) is 
a realization of a PPP with intensity equal to: 

vit) = [ p{t\ s) X{s) ds . (9) 

Jn 

Computing the likelihood p{t\s) requires 
knowledge of the underlying noise properties. We 
will show in section how we determine p{t\ s) 
for our simulated catalogs. 

2.2. Incompleteness 

When a stellar field is observed, some of the 
stars that actually formed within that field cannot 
be detected. They may have evolved into stellar 
remnants (white dwarfs, neutron stars and black 
holes) or been completely destroyed by deflagra¬ 
tion. Some of the remnants may technically be 
observable, but they do not end up in regions of 


® Incompleteness is technically referred to as thinning. 

■^It may be argued that for all modern imaging systems, 
the fundamental observed quantities are fluxes, or better, 
photon counts (or count rates for near infrared detectors). 
Magnitudes are derived quantities and, given the non linear 
relation between fluxes and magnitudes, the noise charac¬ 
teristics are definitely different between them. However, in 
our ideal experiments we assume that we have complete 
knowledge of the noise in magnitude space. Likewise, in 
realistic applications, the noise can be estimated in either 
flux or magnitude space using artificial stars experiments. 
As long as the noise is treated consistently, it shouldn’t 
matter which variable is considered. 
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the CMDj^that can be studied using stellar evolu¬ 
tion models in the traditional sense. We currently 
neglect all evolution beyond the red giant branch 
(RGB) phase in our treatment. The rapid post- 
RGB evolutionary phases could be re-introduced 
into our framework, provided the necessary mod¬ 
els for the later stages of stellar evolution are in¬ 
cluded. 

Other stars might instead not be observable be¬ 
cause they are intrinsically too faint (below the de¬ 
tection threshold) or are at very small angular sep¬ 
aration from brighter neighbors that hamper their 
detection (an effect known as crowding or confu¬ 
sion). Whatever the source of incompleteness, it 
must be accurately modeled and can be included 
in the framework of PPPs. If we indicate with 
0 < a(m, a,z) < 1 the probability of detecting a 
star with a given mass, age and metallicity, given 
a realization ^ = (n, {(m, a, z)i, • • • , (m, a, z)^}), 
the corresponding incomplete realization is ob¬ 
tained by retaining each (m, a, z)i with probability 
a[{in,a, z)i). The incomplete realization is given 
by Ca = (Z, {(m, a, 2 ) 1 , • • • , (to, a, z)i}) with I < n. 
It is possible to show that the incomplete process 
is still a PPP, with intensity: 

Xa{m,a, z) = X{m,a, z) a{m, a, z) (10) 

It is important to notice that incompleteness 
is a property of the models and is a fundamental 
part of the model definition. In section we will 
show a possible way to estimate a{m, a, z). 

2.3. Recap: PPPs for noisy, incomplete 
photometric data and a discrete pa¬ 
rameter space 

Combining everything together, the measure¬ 
ments that one has after observing a patch of the 
sky and performing photometry on the resulting 
images is a realization of an incomplete, noisy PPP 
in a fc-dimensional magnitude space. We are inter¬ 
ested in going from these measurements to a solu¬ 
tion for A, a complete PPP in the 3-dimensional 
space of physical parameters (to, a, z). Solving for 
A means solving for a continuous function over 


®Even if we refer to CMDs and display relevant figures using 
this tool, we emphasize that our method treats magnitudes 
as the real variables, as only individual magnitudes are the 
product of the measurement process, while colors are de¬ 
rived quantities. 


the whole TZ space. This can be accomplished 
by making some simplifying assumptions. First, 
we assume that the IMF is a power-law with a 
single slope 7 , I{m) = ^ oc. nrC. Second, we 
assume that d)(a,z), the combined SFH-MDF, is 
piecewise-constant. In our implementation, the 
constant intervals are evenly spaced; their centers 
form a regular grid. This simplification is equiva¬ 
lent to discretizing the age-metallicity parameter 
space and weighting each grid point by the im¬ 
plied grid-cell’s volume. Solving for A is equiva¬ 
lent to solving for the slope 7 and for the number 
of stars that have formed within each (a, z) cell: 
p( A(to, a,z))= p{ 7 , {rii^ ^} ); for ease of notation, 
ia,z or im,a,z indicate tuples of indices, {ia,iz) and 
{i 7 n,ia,iz), respectively. 

The effective discretization of the parameter 
space implies a substitution of the integral in Eq. 

with a sum over the cells. The final proba¬ 
bility for A, given the incomplete measurements 
77 = (Z,{(Mi,--- ,Mfc)i,--- ,(Mi,--- ,Mfe)i}) is 
thus 


P(7, J h) oc e 





where 



tm,a,z 


}), 

( 11 ) 

( 12 ) 


The indices im,a,z = (*m; iai G) run over the mass, 
age and metallicity cells, and Wi^ is the integral 
of the IMF across the fm-tb cell, which has a 
simple analytic expression in the case of a single 
power law. The index I runs over the observed 
stars. Because the IMF is normalized, we have 
that 


The 1 ^ 01,1 are integrals of incomplete measure¬ 
ment process intensity over the parameter space 
(see Eq. I^and Sect. 2.2). Explicitly, 



lm,a,z 


M 


(13) 
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where 


/ w .,.(0 = 



p{ (Ml, • • • , Mk)i I m, a, z) dm da dz 

(14) 


is the integral of the /-th star’s likelihood function 
over the TZ — cell corresponding to we 


show how f{l) is calculated in Sect. 4.2 


To complete the model specification, we need 
to choose a prior p( 7 , We have already 

specified part of the prior by choosing the range 
covered by the (m,a,z) grid. We discuss the rest 
of the prior specification in Sect. 

Equations (121 and (131 can be modified to add 
further dimensions to the problem, e.g., nuisance 
parameters such as distance or extinction, indi¬ 
cated globally by tt. Adding additional parame¬ 
ters requires defining the incompleteness, stellar 
likelihoods, and intensity prior over the expanded 
parameter set (m, a,z, tt). If we do not attempt 
to solve for the intensity as a function of tt, we 
can marginalize over tt before attempting to solve 
for A(m,a, z); this alters the incompleteness, stel¬ 
lar likelihoods, and intensity prior by a (possibly 
different) constant multiplicative factor at each 
(to, a, z) value. We will show an application in¬ 
volving nuisance parameters in Sect |7.4[ 

We will show in Sect. |4] how the individ¬ 
ual terms of Eq. 11 can be computed. The 


practical solution of the problem of computing 
p{1t {n-ia.A \v) will be given in Sect. ^ where 
we apply a Markov Chain Monte Carlo (MCMC) 
algorithm to generate samples from this distribu¬ 
tion. 


One of the fundamental assumptions underly¬ 
ing our model is that the stars are independently 
independently drawn from the IMF. Such an as¬ 
sumption might not necessarily hold true, accord¬ 
ing to some theories of star formation. In that case 
our method (and all the methods that assume that 
there is an IMF) would fail. The impact of such 
assumption might be more severe in the study of 
young massive clusters, where feedback is partic¬ 
ularly important. However we believe that for the 
study of the Galactic Bulge the independence con¬ 
dition is satisfied. In the Bulge or in other regions 
not actively forming stars, stars that are observed 
within one patch of the sky, have formed possibly 


in different regions, at different times, and have 
undergone dynamical mixing. Therefore, at their 
formation, they were truly independent. If this is 
the case, we can still try and infer the slope of 
an IMF that can be thought as a parent distribu¬ 
tion averaged across the formation history of that 
population. 


3. The model grid 


In order to convert the measured magnitudes 
into physical parameters, i.e., in order to inter¬ 
pret the data in terms of meaningful quantities, 
it is necessary to adopt stellar evolutionary mod¬ 
els. For our examples, we use models computed 


with the Victoria-Regina code (VandenBerg et al. 


20121 , updated with a heavy element mixture 


suited for the stellar populations of the Galac¬ 
tic Bulge (VandenBerg et al. 2014). The trans¬ 
formations from the model physical parameters 
(logL/L©, logTeff, logg) to the observable magni¬ 
tudes are performed using synthetic spectra com¬ 
puted with the MARCS stellar model atmospheres 

2008|). As explained in 


code (Gustafsson et al. 


Sect.|^ we adopt one specific filter system, the one 
for the Wide Field Camera 3 (WFC3) onboard the 
Hubble Space Telescope (HST)] however the scope 
of our method is not limited to one set of stellar 
models or a particular suite of photometric bands. 
Assessment of the systematic uncertainties related 
to the use of different stellar models, as well as the 
exploration of the best possible combination of fil¬ 
ters for deriving IMF, SFH and MDF of resolved 
stellar populations is beyond the scope of this pa¬ 
per. 


For this work, we use a grid with a 2% spacing 
in mass, 2.5% in age, and O.I dex in [Fe/H]. The 
mass and age steps correspond to a constant spac¬ 
ing in the logarithm of mass and age, respectively. 
This type of spacing translates into a conveniently 
more constant spacing between stellar models in 
the CMD. However, the grid spacing does not 
have to be regular for our technique to be applica¬ 
ble. Given that this method was developed to deal 
with multi-band photometric data for the Galac¬ 
tic bulge, we choose a range of parameters that is 
suitable for the bulge stellar population; again, the 
grid range can be changed and adapted to different 
problems. Specifically, we have to S [0.2,1.5] Mq, 
a e [7., 14.7] Gyr, z € [-2.0, -bO.4] dex. 
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The choice of the grid size and grid resolution 
has a practical impact on the method and is a 
matter of compromise between computing time 
and the ability to resolve rapid changes in $(a, z). 
The number of unknowns in our problem is al¬ 
most equal to the number of age-metallicity cells. 
Because we choose to parametrize the IMF as a 
single power law, one additional parameter is the 
slope of the IMF, 7 . More complex models for 
the IMF could in principle be chosen in other ap¬ 
plications, if there is reason to think that -or to 
explore whether- the data may help constrain an 
IMF model with increased complexity. 

The choice of the mass, age, and metallicity 
grid resolution must also be guided by the qual¬ 
ity of the data at hand. One should be able to 
resolve the individual likelihoods without compu¬ 
tationally expensive excessive oversampling. Un¬ 
fortunately, stars in different evolutionary phases 
have very different likelihoods, in terms of how dif¬ 
fuse each likelihood is in the physical parameter 
space (m,a, z). For example, if we can resolve a 
turnoff star’s age likelihood function, we are over- 
sampling the main sequence stars’ age likelihood 
functions. We can use the typical photometric er¬ 
rors in the turnoff region to estimate the precision 
attainable in age determination and then use a 
fraction of that as the age-step of the grid. The 
use of metallicity-sensitive photometric bands can 
increase the resolving power of our observations 
for [Fe/H]. The WFC3/UVIS F390IU band, for 
example, covers a spectral region that, in dwarf 
stars, contains strong metal lines. 

4. Noise, incompleteness and likelihoods. 

We test the effectiveness of our method by simu¬ 
lating data and comparing the input and recovered 
intensity functions. The simulated quantities are 
stellar magnitudes. We use Gaussian noise for the 
simulated magnitude measurements and assume 
that the only errors involved are random errors. 
We further assume that the noise in multiple pho¬ 
tometric bands for the same star is uncorrelated. 
For catalog generation and incompleteness evalu¬ 
ation, we also assume that the detections in each 
band are independent of each other, and that a 
star is considered detected when it is detected in¬ 
dependently in all bands. 

To estimate realistic noise levels for a typical 


stellar population study, we use data from Brown 


et al. (2009 2010). These papers present obser¬ 


vations of four stellar fields towards the Galactic 
bulge. The data from the full observing program 
will be analyzed in a forthcoming paper; for test¬ 
ing our method, we use the subset of the data that 
was taken in the OGLE29 low-extinction window. 
We also limit ourselves to 3 of the 5 bands, namely 
F390W, F555W, and F814W, respectively approx¬ 
imating Washington G, Johnson V, and Johnson 
/, which is sufficient to demonstrate our method. 
Because our simulations do not include distance 
and extinction as parameters, we convert the ap¬ 
parent magnitudes into absolute magnitudes by 
subtracting the average bulge distance modulus 
and typical extinction in each of the 3 bands. 

For each band, we use the real data to esti¬ 
mate a relation between the observed magnitudes 
and photometric uncertainties. We first bin the 
observed stars in magnitude bins and take the av¬ 
erage photometric uncertainty in each bin. The 
uncertainties come from PSF fitting photometry 
using DAOPHOT (Stetson 1987). We then fit a third- 
order polynomial to the logarithm of the average 
uncertainty as a function of bin magnitude. The 
fit coefficients are then used to compute the typi¬ 
cal photometric error as a function of magnitude 
(see Fig. [^. 

We also adopt an incompleteness vs. mag¬ 
nitude relation. This relation is obtained by 
imposing that objects with magnitude errors of 
0.05,0.1,0.5, and 1 mag are detected in 100%, 
95%, 50%, and 0% of the cases, respectively, and 
linearly interpolating the incompleteness between 
these magnitude values. These incompleteness 
curves are shown in Fig. [T| As in the data of 
Brown et al.| ( |200^[Mot , the F555W and F814W 


bands have lower uncertainties and incomplete¬ 
nesses at fixed magnitude than the F390W band. 
This means that the overall incompleteness of the 
simulated observations is generally dominated by 
F390W. There could be situations in which com¬ 
pleteness does not only depend on the stellar mag¬ 
nitudes, but also on the position across the field of 
view. This is specially true for dense stellar clus¬ 
ters, given the gradient in the crowding properties 


between the center and outskirts (see e.g. Gennaro 


et al. 2011). However the present work has been 


developed to deal with HST observations of the 
Galactic Bulge (Brown et al. 2009 2010), with 




























a small field of view (~ 4 arcmin) across which 
crowding is very homogeneous. We will limit our¬ 
selves to this simple situation where incomplete¬ 
ness does not depend on position. 

4.1. Computing a(TO, a, z) 

The incompleteness curves of Fig. represent 
the probability of detecting a star given its mea¬ 
sured magnitude. The function a(jn,a,z) repre¬ 
sents the probability of detecting a model star, 
for which we know the intrinsic, or error-free, 
magnitudes (M 390 , M 555 , M 814 At a single 
value of (m, a, z), the model incompleteness is an 
average of the measured-magnitude incomplete¬ 
ness over the measured-magnitude pdf. Since 
we require a star to be detected in all bands, 
the incompleteness of a vector of measurements 
(Mago, M 555 , M 814 )°’^® is the product of the incom¬ 
pleteness in each band. 

To compute a for each model grid cell cen¬ 
tered on {m,a,z), we simulate the attempted 
measurement oi j = 1 , • • • , 1000 stars per grid 
cell. The 1000 (m, a, z)™*’’ values are randomly 
uniformly extracted within the cell. For each 
of them, (M390, M555, M814)™*"' is computed us¬ 
ing the library of stellar models. Corresponding 
(M390, M555, M814)°'^® are extracted from Gaus- 
sians centered on (M390, M555, M8i4)™*‘' with the 
appropriate tr’s. Finally, using the curves in Fig.[T] 
we decide whether that star would have been ob¬ 
served or not. This is done by comparing, for each 
band, the value of the detection probability with 
a uniform random number between 0 and 1. If a 
star is not detected in one band, then it is consid¬ 
ered not detected at all. The number of recovered 
stars divided by 1000 is the approximate value of 
a(m, a, z) averaged over the model cell. 

There are several simplifying assumptions in 
our treatment. In general, for real observations, 
the explicit form of the noise is not known (even 
though it is often assumed to be Gaussian). How¬ 
ever, in the case of real observations, where real 
images of stellar fields are used, the process de¬ 
scribed above can be reproduced using artificial 
star tests. The latter consist of introducing into 
the images under study stars of know magnitudes 
(directly related to known (m, a, z) through stel¬ 
lar models) and then trying to measure the output 
magnitudes. By using the exact same algorithmic 


sequence (including cuts, a posteriori selections, 
detection criteria and so on) as used for building 
the real catalog of observations, we can prepare 
a list of input vs. output (detected and unde¬ 
tected) stars. For the computation of a{m,a,z) 
it is still possible to generate (m, a, z)™*"' and the 
corresponding (M3go, M555, M814)™*'' from stellar 
models. The list of input stars can be searched for 
stars with these magnitudes and the correspond¬ 
ing output can be checked to see whether it cor¬ 
responds to a valid detection. Finally, the ratio of 
detection-to-input stars can be used as a(m, a, z). 

The above description deals with incomplete¬ 
ness due to the measurement process and crowd¬ 
ing. These phenomena mostly affect models to¬ 
wards the faint (low-mass) end of the parameters 
grid. However, for the treatment of PPPs, we also 
need to consider stars that contribute to A but are 
no longer observable because they have evolved 
away from the main sequence. These stars are in¬ 
stead towards the high-mass end of the model grid. 
To account for these lost stars, it suffices to check 
whether the (m, a, z)™**^ points that are generated 
per grid cell in the previous step correspond to 
“alive” stars. If they do not then they are simply 
considered as non-detected models and contribute 
to lowering a for that cell. 

For real observations there might be cases 
where further a posteriori cuts are imposed to 
photometric catalogs. For example, stars near 
the detection threshold are often removed because 
their measurement uncertainties are difficult to 
estimate. 

In our simplified scheme, we assume we know 
everything about the noise model and we also men¬ 
tion that, in realistic cases, artificial star tests can 
be used to explore the noise characteristics and 
derive the incompleteness function. However, even 
artificial star experiments may not capture all pos¬ 
sible sources of errors and, in general, one might 
be wary of trusting the noisiest detections. If this 
is the case, any hard magnitude cut imposed to 
a catalog must be included in the computation 
of aim, a, z). The only difference with the pre¬ 
viously illustrated scheme is that the simulated 
values (M390, M555, M814)j'^® must also fall above 
the imposed thresholds in order to be accounted 
for as detections. We will show the impact of such 
cuts on the recovery of A in Sect. 
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Fig. 1 .— Realistic values of the typical error and completeness derived from data by Brown et al. (2009 


20101 , for the three photometric bands considered in this work. The x-axis indicate the measured magnitudes. 


Solid lines: error curves; dashed lines: completeness (detection probability); dotted vertical lines: 50 % 
completeness limit for that band. 


Figure shows a{m,a,z) for our grid of mod¬ 
els assuming the standard curves of Fig. in¬ 
completeness is color-coded from dark blue ( 
a{m,a,z) = 0 ) to red ( a{m,a,z) = 1 ) . Each 
model point corresponds to the center of a (m, a, z) 
cell. The CMDs in the top row show the values of 
a{m, a, z) when no cut is applied, while in the mid¬ 
dle row, a cut at Mg,i4 = 9 mag was applied, hence 
the sharp drop in incompleteness. To illustrate the 
effects of stellar evolution on incompleteness we 
separate the contribution to a from stars that 
have evolved, shown in the bottom row; here it is 
possible to notice that for the most massive model 
cells (in the red giant phase) a fraction of the cell 
volume corresponds to stars that are no longer 
observable, even when the central values (m, a, z) 
for that cell correspond to still-alive stars. For our 
method, the total a{m, a, z) is the product of the 
incompleteness in the first (or second) row with 
the incompleteness in the third row. 


4.2. Individual stellar likelihoods 


An individual likelihood is the probability of de¬ 
tecting a star at (M390,M555,M8i4)°'^® given the 
model parameters (m,a, z). The discretization 
of Eq. ( 13 ) implies that the likelihoods, fi^ ^ ^{l) 
need to oe integrated across the individual cells 
of the parameters grid. To compute f{l), we use 
the following method. Eor each parameter grid 


cell we generate 500 values of (m, a, z)j within the 
cell. We compute the corresponding magnitudes 
(M390, M555, M814)™°‘^ using the stellar models li¬ 
brary. We then evaluate the photometric a’s at 
those magnitudes. The magnitudes and cr’s com¬ 
pletely characterize the Gaussian likelihoods; for 
each star i in a simulated catalog, with magnitudes 
(M390, M555, M814)°’^® we compute the value of 


p((M39o, M555, M 814 )°*’'* I irn,a,z)j) = 



'\/27r(Tfj[t,4 


X exp(-xV2) 


( 15 ) 


with 


X 


2 


y / 

m= V 

390 , 555,814 


2 


( 16 ) 


We then average the likelihood over the 500 j- 
values. This corresponds to a marginalization of 
the likelihoods across the model grid cell. 

The necessity of averaging the likelihoods over 
a model cell, instead of just taking the likelihood 
value for, e.g., the cell center, is apparent when 
considering evolutionary phases in which a star’s 
position on the CMD rapidly changes. While on 
the main sequence the changes are very slow, once 
a star reaches the turnoff region it may move a 
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Fig. 2.— CMDs of our model grids, color coded by incompleteness values. Red corresponds to 100% complete 
models, blue to 0%. The first two rows show the effect of incompleteness due to missing detections when 
no a posteriori cut is applied to the photometric catalogs (top row) or when a cut at Msi 4 = 9 mag is 
applied (center row). The bottom row shows the contribution to incompleteness from stellar evolution. Red 
corresponds to model cells that are not evolved away from the main sequence; bluer colors correspond to 
cells where some fraction of the models is not observable, having evolved away from the main sequence. 
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Fig. 3.— Illustration of our likelihood calculation method based on averages over the cells. The blue 
points represent 500 model realizations within cells centered on M = 0.59Mq (left), andO.92M0, (right). 
Both model cells are centered on a = 12.05 Gyr and [Fe/H] = —0.3 dex. The black points are simulated 
measurements from one of our artihcial catalogs (Blg, see Sect.|^. The pink points are data that are within 
5cr from at least one of the 500 realizations (where a is the uncertainty associated to each model realization, 
according to the curves of Fig. [^. The shade goes from high (light pink) to low likelihood (dark pink). The 
error symbols correspond to the models with parameters in the center of the grid cells, and their size is equal 
to the typical a of that cell. 


substantial length (in both magnitude and color) 
within a very short time. Because these fast 
phases happen at different times for different stel¬ 
lar masses and metallicities, it would be compu¬ 
tationally very intensive to define the models on 
a grid that is fine enough to resolve the fastest 
stellar evolutionary phases for all (m, a, z) com¬ 
binations. However, the average across a cell al¬ 
lows for a proper treatment of these phases even 
within the limits of the grid resolution. Figure 
shows an example of such a situation for one of our 
simulated catalogs (black points), focusing on two 
different regions of the CMD: the main sequence 
on the left and the turnoff on the right. In light- 
blue we show the 500 models realizations for two 
different cells, centered on M = 0.59 and 0.92Mq 
respectively. Both cells are centered on a = 12.05 
Gyr and [Fe/H] = —0.3 dex. The width of the cells 
is 2% in mass, 2.5% in age and 0.1 dex in metal- 
licity. The light-blue area is the transformation of 
the physical parameter cube (m±(5m, azLSa, z±5z) 
into the CMD. 

The models at the cell centers are identified by 
the error symbols, with error bar size equal to the 


average a on the cell. The symbols in pink are 
simulated catalog stars that are within 5 (t of at 
least one of the 500 realizations (where cr is the 
uncertainty associated to each model realization, 
according to the curves of Fig. [^. The shades 
of pink correspond to high-to-low (light-to-dark) 
likelihood values on a logarithmic scale. The need 
to compute likelihoods over entire model cells is 
particularly clear in the right panel. If we had 
considered only the model corresponding to the 
grid cell center, many data points would have been 
too far (in units of cr) to have any likelihood at that 
cell even though their physical properties are very 
close to the physical properties of the cell’s center. 
Averaging the likelihood over the entire model cell 
helps prevent this problem. 

Figure shows how the stellar parameters are 
constrained differently in different stellar phases. 
The top row shows an RGB star with a mildly 
constrained age and very well constrained metal- 
licity. The middle row shows a turnoff star for 
which the likelihood is very narrow in all direc¬ 
tions. Finally, in the bottom row we show a main- 
sequence star, for which mass and metallicity are 
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constrained moderately well, but for which age is 
almost completely unconstrained, at least within 
the range of ages adopted here. 


5. Solving for the intensity function 


As detailed in Sect. |2.3[ the problem of solv¬ 
ing for A can be regarded as solving for the IMF 
slope ( 7 ) and the number of stars formed in each 
age-metallicity bin, Equation (11) de¬ 


scribes the posterior pdf from which we wish to 
sample in order to obtain a solution. The gen¬ 
eral shape of star-formation histories in the age- 
metallicity plane is expected a-priori to be sparse. 
This means that only a limited number of cells 
have occupancy greater than 0. In actual 

stellar fields, for example, the existence of an age- 
metallicity relation will greatly reduce the num¬ 
ber of active model cells. There are existing tech¬ 
niques, developed for image reconstruction, that 
allow one to deal with situations in which only a 
few pixels of the image contain the desired infor¬ 
mation while the rest of the pixels constitute a 
noisy background. These techniques allow one to 
suppress the background and sharpen the signal 
in the pixels where it is present. To improve our 
solution, we adopt an approach similar to that de¬ 


scribed in Lingenfelter et al. (2009) and designed 
for Poisson-distributed data. This approach is 
equivalent to imposing a Lomax, or Pareto Type 
II, distribution as the prior (see Eq. II) on the 
pixel intensities rii^ ^. The product of all of the 
cells’ priors is proportional to 


n( 




1 


(17) 


In the language of regularization, /? sets the thresh¬ 
old, or minimum value that is considered to be 
an actual, rather than spurious, signal, and S sets 
the threshold’s sharpness, or strength with which 
values below the threshold are drawn towards 
zero. After experimenting, we adopted 5 = 2 and 
/3 = 0.4, meaning a soft threshold, and meaning 
that we regularize (put to 0 ) only (a, z) model cells 
with very low occupancy. At this value of /3, the 
Lomax distribution is improper, i.e. integrates to 
infinity. Because the likelihood function of a PPP 
decays exponentially at high values of A, the pos¬ 
terior pdf is still proper. 

For the IMF slope 7 , we assume a uniform 


prior bounded between —3 and -|-3. This range 
includes all of the commonly assumed high-mass 
IMF slopes. 

Given the high dimensionality of this prob¬ 
lem [D = 776 with the grid used here), sampling 
the whole space of ( 7 , ^}) can be inefficient, 

implying slow convergence when using traditional 
Markov Chain Monte Carlo sampling methods. In 
order to accelerate convergence, we first estimate 
the maximum-a-posteriori (MAP) solution for 
(A, {rii^^}) using an Expectation-Maximization 
algorithm (see Appendix 0 . We then start from 
the MAP estimate and use Metropolis-Hastings 
(M-H) Markov Chain Monte Carlo (MCMC) to 
generate samples from the posterior pdf. M-H 
MCMC performs adequately and has the addi¬ 
tional merit of being simple and intuitive. 

6. Simulated catalogs 

To test our method, we simulate catalogs 
with a range of star-formation histories and age- 
metallicity relations. We explore the impact of 
increased noise levels in the data and the impact 
of different magnitude cuts applied to the cat¬ 
alogs. We also show how results differ when the 
metallicity-sensitive band 7^39011^ is removed from 
the analysis. 

For each catalog, we define a shape for A(m, a, z) 
X{m) $(a, z). We also specify the true number of 
expected stars (A^true), equivalent to specifying the 
integral of A. The actual number of stars formed 
(A^born) is extracted from a Poisson distribution 
with expected value equal to A^true- Given the 
IMF single power-law parametrization, we need 
only specify the slope parameter; we set 7 = — 2 
for all the simulated catalogs. The function <i)(a, z) 
represents the age-metallicity relation for all the 
stars that are formed. We generate catalogs with 
isolated narrow bursts of star formation, as well as 
more extended star-format ion histories, with large 
changes of stellar metallicity with age. Table 
summarizes the characteristics of each catalog. 

Each star’s (m, a, z) values are extracted from 
I(m) and 4)(a, z). Given the finite lifetime of 
stars, not all of the extracted values correspond to 
observable stars. Stars that are too massive for the 
extracted age and metallicity are still considered 
for the budget of formed stars but they obviously 
do not end up in the observable catalog. If instead 
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Fig. 4.— Examples of likelihoods for stars in different evolutionary phases: RGB (top row), turnoff {middle 
row), and main sequence {bottom row). The left panels show the stellar location in the CMD. The other 
3 panels show the likelihoods as functions of different pairs of variables, marginalized over the third. The 
orange circles indicate the observed magnitudes {left column) and the true mass, age, and metallicty values 
{columns 2-4). 


the extracted age does not exceed the maximum 
lifetime at given m and z, stellar models are used 
to compute the intrinsic stellar magnitudes. We 
then use our noise model to assign each star a set 
of measurable magnitudes, and include or exclude 
the star from the final catalog based on the in¬ 
completeness curves in Fig. [T] For some of the 
catalogs, we use a version of the error curves in 
Fig.g where the whole curve is multiplied by a 
factor (/err) equal to 2.5, in order to explore the 
effects of reduced data quality. The completeness 
curves are also recomputed accordingly. Figure 
shows the extracted parameters for catalog Blg. 
The histograms represent all of the formed stars 
(light blue) with the subset of observed stars su¬ 
perimposed (dark blue). The contour plot shows 
the form of <i)(a, z) chosen for this catalog. The 
three panels in the bottom row of Fig. show the 
CMDs for the same catalog. Analogous figures for 
the other catalogs are shown in Appendix [B| 


7. Results 


We use the outputs of the MCMC runs, Sect.[^ 
to define our estimates of the parameters and 
their uncertainties. The results are summarized 
in Figs. |6 - 11 In all the panels, different shades of 
blue represent the recovered solutions, while the 
original values are in orange. 

For the 2D star-formation history (left, bot¬ 
tom), we show the MAP solution as a way to 
represent the best combination of {ui^ values. 
Defining the uncertainty interval for the recov¬ 
ered solution in the age-metallicty plane is not 
straightforward. While we have samples for the 
ill it is the global solution for the whole (a,z) 
plane, and its uncertainty, that is of interest. Ad¬ 
jacent cells in the (a,z) grid are strongly coupled 
by the fact that most of the individual stellar like¬ 
lihoods are spread over large age-metallicity re¬ 
gions. This means that the pdf for the 
has a complex covariance structure, which is hard 
to represent and convey in a meaningful way. In¬ 
stead, we show the uncertainties of different col- 
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Table 1 

Main characteristics of the simulated catalogs 


Name 

7" 

4>(a,z) 

true 

^born 

Nobs 

/err " 

^ d 

CTAMR 

Bug 

-2 

Bulge-like SFH and MDF 

25000 

25256 

9921 

1 

0.05 

Blg:No390 

-2 

Same as Blg, no F390W data 

25000 

24910 

18480 

1 

0.05 

BlG:F2p5 

-2 

Same as Blg, larger error 

25000 

25123 

6264 

2.5 

0.05 

BlgiLrg 

-2 

Same as Blg, more stars 

100000 

100138 

39246 

1 

0.05 

Burst3 

-2 

Three equal-intensity bursts 

25000 

25002 

10256 

1 

0. 

Exp 

-2 

Exponential decay 

25000 

25237 

12298 

1 

0. 

Const 

-2 

Constant star formation 

25000 

25404 

11624 

1 

0. 


Note.— The BLG catalogs have SFH and MDF that emulate observations of the Galactic Bulge. 
The Bursts catalog has burst at (a, z) = (12.5, —1.0); (8.5, —0.3); (8.5, —1.0). In both the Exp and 
Gonst catalog there is a linear age-metallicity relation with a € [10,12] Gyr and ^ € [—1.5, —0.5] 
(see also Figs. 13 and 141. The exponential decay is such that within the overall time interval, the 
star formation decays by 3 e-folds. 


““IMF slope 

’^Combination of star-formation history and age-metallicity relation 
‘^Multiplication factor for the error curves with respect to the standard curve of Fig. 
‘’Dispersion about the mean age-metallicty relation 


lapsed 1-dimensional versions of We plot 

the marginal age (left, top) and marginal metal- 
licity (center, bottom) distributions, as well as the 
cumulative star-formation history (CSFH) (right, 
bottom) and the mean age-metallicity relation 
(AMR) (right, top). 

For the first three of these collapsed distribu¬ 
tions, we use the geometric medians as our best 
estimates. The geometric median is a commonly 
used central tendency indicator; it is the mul¬ 
tidimensional generalization of the common me¬ 
dian. In this case, we want to identify a me¬ 
dian marginal SFH, or MDF, or CSFH, and each 
grid point where such distributions are computed 
constitutes a dimension. We treat the marginal 
pdf, for one MCMC iteration, as a point in an 
TOgrid-dimensional space, where Wgrid is the num¬ 
ber of grid points along the age or metallicity di¬ 
rection. When considering all the MCMC itera¬ 
tions, the geometric median is the point in the 
TOgrid-dimensional space that minimizes the sum 
of the (Euclidean) distances to all the other points. 
To quantify the uncertainty on the marginal dis¬ 
tributions, we consider each bin individually. We 


take the difference between the bin value of the 
geometric median and the bin value of the distri¬ 
butions for each MCMC iteration. We then sort 
the differences and take the 16% (84%) quantiles; 
these quantiles are analogous to the typical 1 — cr 
interval in case of a Gaussian distribution. 

We define the mean AMR as the mean value of 
the metallicity of stars formed within one age grid 
cell. The AMR can have an intrinsic dispersion (as 
in the case of all the BLG-type catalogs), meaning 
that at fixed age, a star can be formed with a range 
of metallicities, as in the case of non-instantaneous 
mixing of the ISM. Given the possible intrinsic 
dispersion, in the summary plots we only show 
the comparison between the simulated AMR and 
the MAP solution. To avoid confusion, we do not 
show the additional dispersion introduced in the 
solution by the fact that each MCMC sample can 
have a slightly different mean AMR. It is worth 
noting that in most cases the AMR solutions in 
the top-right panels depart significantly from the 
input AMR (orange points). At a closer look, how¬ 
ever, it is clear that the departure is limited to 
metallicities where the 2D solutions are generally 
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Fig. 5. — Generated parameters (top) and simulated photometry (bottom) for the Blg catalog. The contour 
plot shows the form of $(a, z). The three histograms show the mass, age and metallicity distribution for the 
whole simulated population (light blue) and the stars that make it to the simulated, observed catalog (dark 
blue). The latter are shown in three different CMD combinations in the bottom panels. 


negligible. For these metallicities the input AMR 
is not defined (no stars are formed in the input 
model), but the output mean AMR can still be 
computed. When looking simultaneously at the 
AMR and 2D solution, it is clear that the recov¬ 


ered AMR could be, in fact, truncated to metallic¬ 
ities where the 2D solution is significant. However 
we show the full derived AMR for the sake of com¬ 
pleteness, 

Finally, the IMF slope panel (center, top) shows 
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Fig. 6. — Recovered properties for catalog Blg. Left, bottom: the maximum-a-posteriori solution for 
{riia ^}. Left, top: star-formation history. Center, bottom: metallicity distribution function. Center, top: 
distribution of the IMF slope values from the MCMC runs. Right, bottom: cumulative star-format ion 
history. Right, top: mean age-metallicity relation and its dispersion. In all panels the simulated values are 
in orange, the recovered ones in shades of blue. 


the input value of (orange) and posterior pdf for 
the IMF slope 7 . 

Figurej^shows the results for catalog Blg. The 
distribution of the recovered IMF slope is cen¬ 
tered on the true value, and the marginal SFH 
and MDF match the input. Some deviation is ob¬ 
served at old ages for the SFH, CSFH, and AMR, 
most likely because few stars formed at these ages 
and even fewer were observed. Moreover, the stars 
are formed at low metallicity; the isochrone colors 
become increasingly degenerate at low metallici- 
ties, leading to individual stellar likelihoods that 
are less well-defined, and looser AMR constraints. 

7.1. Catalog cuts and photometric errors 

In the case of Fig. we have used all of the 
observed stars of catalog Blg to reconstruct the 
intensity function. However, as anticipated in 
Sect. |4.1[ there might be cases in which it is de¬ 
sirable to adopt conservative cuts, avoiding stars 
at the faint limit. As illustrated in Sect. |4.l| these 
hard cuts can be accounted for in the incomplete¬ 
ness function since they are simply another aspect 
of the detection process. We show in Fig. [^the ef¬ 


fect of cuts at Msi 4 = 7 (top) and 9 (center) mag, 
respectively. The = 9 mag cut corresponds 
approximately to a cut at the 50% completeness 
limit (see Fig. [^. 

Such cuts affect different aspects of the solution 
in different ways. Generally speaking, the shape 
of the SFH, MDF, and AMR are not significantly 
changed. This is to be expected, since both of 
the applied cuts leave the turnoff and RGB intact. 
The turnoff carries most of the age information, 
while the RGB is a good metallicity indicator. The 
cuts have some effect on the details of the SFH and 
MDF, especially in poorly-populated parts of the 
(a, z) plane and where there are few, if any, turnoff 
or RGB stars. The largest difference is in the IMF 
estimates. The reduced mass range implies that 
the IMF slope cannot be recovered with the same 
accuracy and precision. This in turn affects the 
overall normalization of the SFH, as evident in 
the CSFH plot in the = 7 mag cut case. It 
is encouraging that the Msu = 9 mag cut case 
still looks very good, in all respects. Although the 
precision of the 7 recovery is lower than in the full 
catalog case, there is no obvious bias (in contrast 
to the Msi 4 = 7 mag cut case). 
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In an ideal situation, where one has perfect 
knowledge of the completeness and photometric 
errors at all magnitudes, it would be best to ana¬ 
lyze the full catalog. With real data, one must be 
careful using stars near the faint limit, because the 
effects of any biases in the characterization of the 
completeness or photometric uncertainties will be 
amplified. Fortunately, a conservative cut to the 
catalog at the 50% completeness limit, with proper 
inclusion of this cut in the selection function, can 
yield results that compare well to those from the 
full catalog. For the rest of the examples, we will 
be mostly using catalogs with a M 814 = 9 cut. 

The bottom panel of Fig.j^shows the results for 
catalog Blg:F2p5, simulated with the same IMF, 
SFH, and MDF as catalog Blg, but with photo¬ 
metric errors increased by a factor of 2.5; we limit 
the comparison to the case with M 814 = 9 mag 
cut. The overall recovery in this case is similar 
to the Blg case. In particular, the IMF recov¬ 
ery looks very similar. This is because individ¬ 
ual stellar masses are still very well recovered, at 
M 814 < 9, because our mass grid spacing (2%) 
is coarse enough that the Blg and Blg:F2p5 
cases are almost indistinguishable. Similarly, the 
marginal (and cumulative) SFH and MDF are also 
well recovered. However, the details of the AMR 
start to get worse for the age bins with fewer stars; 
with the larger uncertainties, the information con¬ 
tent of the few stars that were generated in those 
bins is no longer sufficiently constraining. This ex¬ 
ample demonstrates the need, in real cases, to un¬ 
derstand the limitations of the available dataset. 

7.2. Impact of Nbands and Nobs 

Catalog Blg:No390 has been simulated using 
the same IMF, SFH, and AMR as catalog Blg, 
but in this case we have changed the mock ob¬ 
serving strategy to only consider 2 filters, omit¬ 
ting F390W. This choice enables deeper catalogs 
(see the number of stars in Table at the cost 
of having weaker metallicity constraints. In our 
simulations, F390W is the shallowest band and 
has the worst signal-to-noise ratio. This is a com¬ 
mon choice when designing an observing strategy: 
more filters and a shallower catalog vs. deeper 
observations and less chromatic information. The 
top panel of Figure [^demonstrates that the deep 
observations guarantee a good result for the IMF 
slope, and hence the normalization of the SFH. 


The mean AMR, SFR, and MDF match the input 
to within the uncertainties. However, the MDF 
itself is highly uncertain in this observing config¬ 
uration, as expected. 

The bottom panel of the same figure shows 
the results for catalog Blg:Lrg, whose $(«, z) is 
equal to 4 times that of Blg; this corresponds 
to observing the same population over an area 
4 times larger. The recovery of all quantities is 
nearly optimal. 

7.3. Differentiating similar histories 

In some real cases it might be very important 
to be able to distinguish between apparently sim¬ 
ilar star formation scenarios. Catalogs Exp and 
Const have been designed to test the ability of 
our method to differentiate such cases. The cat¬ 
alogs have the same IMF, AMR, photometric er¬ 
rors, and selection criteria, but the star-format ion 
rate is exponentially decaying in one case and con¬ 
stant in the second. The total number of stars 
formed and the duration of star formation are also 
the same. The scenarios’ CMDs, which are shown 
in Figs. [Island |l4l appear quite similar. 

The results for the recovery are shown in Fig.[^ 
for the cases in which both catalogs are cut at 
Af 8 i 4 = 9 mag. The recovery is again largely 
successful, with a clear distinction between the 
two scenarios. The differences are clearest in the 
CSFH plot. 

7.4. Nuisance parameters: the effect of a 
distance distribution 

As was explained at the end of Sect. |2.3[ our 
model can incorporate nuisance parameters (NPs) 
at the cost of increased computational time. We 
have tested the effects of introducing NPs by simu¬ 
lating catalogs with the same I{m) and $(a, z) as 
those in Tab.[^ but with a distribution in distance 
modulus (DM). In particular, we will show results 
for the Burst3 catalog case, where three isolated 
peaks of star formation were simulated. This cat¬ 
alog provides a good visualization for the effects 
of a distance spread. With respect to the original 
catalog, the one including a distance distribution 
has been spread using a Gaussian distribution in 
DM with expected value 0 and ctdm = 0.25 mag. 
This dispersion is similar to the spread in DM of 
stars along Galactic bulge sightlines. 
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Fig. 7.— Recovered properties for catalog Blg, with cut at = 7 mag (top), M 814 = 9 (center), and 
catalog Blg:F2p5, cut at M 814 = 9 mag (bottom). See Fig. |^for a description of the individual panels. 
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Fig. 8. — Recovered properties for catalog Blg:No390 (top) and catalog Blg:Lrg (bottom), both cut at 
Msia = 9 mag. See Fig. for a description of the individual panels. 


As outlined at the end of Sect. 2.3 when dealing 
with nuisance parameters both the incompleteness 
and likelihood function have to be computed as a 
function of the nuisance parameters too. A prior 
must then be specified in order to marginalize over 
them. 


In Fig. |10[ we show results where distinct pri¬ 
ors are assumed for the distance distribution (one 
correct and two incorrect). Specihcally, the re¬ 
sults in the three panels are obtained under the 
following assumptions: a Gaussian prior with ex¬ 
pected value 0 and cdm = 0.25 mag (second from 
the top; the correct prior), a uniform prior with 0 
mean and 0.25 mag width (second from the bot¬ 
tom; an incorrect prior that is still a reasonable 


approximation), and a single distance centered on 
0 (bottom; an incorrect prior that is a less rea¬ 
sonable approximation). For comparison, we also 
show the recovery in the case where all the stars 
are simulated at the same distance (top). 

With a spread in distance, the degradation in 
the recovered IMF, AMR, SFH, and MDF is quite 
obvious. Even in the case of the correct prior, 
the fact that the real distances of individual stars 
are known only with some probability makes both 
the individual ages and metallicities more uncer¬ 
tain. Essentially, all of the distributions have been 
convolved with the corresponding distance uncer¬ 
tainty. However, at least in the case with a cor¬ 
rect distance prior (second from the top), the re- 
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Fig. 9. — Recovered properties for catalog Exp (top) and catalog Const (bottom), both cut at = 9 
mag. See Fig. for a description of the individual panels. 


covered values match the input to within uncer¬ 
tainties. In contrast, incorrect assumptions on the 
distance prior further increase the discrepancies 
between the truth and the result. In particular, 
assuming that all the stars are at the same dis¬ 
tance makes all the answers for AMR, SFH, and 
MDF completely wrong while the IMF slope is less 
affected. 

These examples constitute a serious warn¬ 
ing against oversimplifying assumptions when it 
comes to priors. In the case of the Galactic bulge, 
it is certainly true that the DM dispersion is non- 
negligible, and so particular care should be taken 
regarding the assumed DM distribution. 


7.5. Applying spectroscopic constraints 

Often, in the study of resolved stellar pop¬ 
ulations, there are additional constraints that 
can be useful when solving for their IMF, SFH, 
and MDF. Spectroscopic constraints, usually from 
RGB stars, can improve the solution if they are 
appropriately handled. As is the case for pho¬ 
tometric incompleteness, it is very important to 
have a good knowledge of the selection function 
that is used to build the sample of targets for the 
spectroscopic observations. 

A catalog of spectroscopic measurements and 
errors can be added to the catalog of photometric 
measurements to further constrain the intensity 
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Fig. 10. — Recovered properties for catalog BurstS. Top: the case without distance spread. Second from 
the top: the correct distance distribution is used as distance prior in the recovery. Second from the bottom: 
a uniform distance prior between —1 and +1 (t is used instead of the correct (Gaussian) one. Bottom: all the 
stars are considered to be a — priori at the same distance, equal to the mean of the correct Gaussian prior. 
In all cases we adopted a catalog cut at Msi 4 = 9 mag. See Fig. for a description of the individual panels. 
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function (A). Technically, the spectroscopic data 
constitute a new, related PPP, in the spectroscopic 
measurement space, with the same (up to nor¬ 
malization) underlying true intensity as the one 
underlying the photometric measurements PPP. 
The probability of the combined spectroscopic and 
photometric PPPs is the product of the probabil¬ 
ities of each one: 


p{{vspec} J {??phot} I A) — 

P({^spec} I A, A^spec) ^ P({^phot} | A), 


where A^spec is the (possibly different) normaliza¬ 
tion of the spectroscopic PPP. 


The selection function for the spectroscopic 
sample can be regarded in the same way as the 
incompleteness function of photometric data and 
computed with the same methods described in 
Sect. |4.1[ Once this function is evaluated, we can 


write an equivalent of Eq. (13) for the spectro¬ 
scopic PPP. 


We simulate the process of building a spectro¬ 
scopic sample and obtaining corresponding mea¬ 
surements by first generating a photometric cata¬ 
log with the same true properties as the Blg cat¬ 
alog, but a larger total number of stars - as noted 
above, the normalizations of the spectroscopic and 
photometric PPPs do not have to be same. We 
then select RGB stars with 2 < Msi 4 < 3 mag. 
This is analogous to selecting targets from a shal¬ 
low wide-area survey centered on a similar posi¬ 
tion as the field where deep imaging is available 
for IMF, SFH, and MDF reconstruction. We then 
assign to each spectroscopic target an [Fe/H] error, 
extracted from a gamma distribution with shape 
parameter 50 and scale parameter 0.001, thus ob¬ 
taining a mean error of 0.05 dex, and a standard 
deviation of the errors 0.007 dex. Finally, for each 
star we extract the measured [Fe/H] from a Gaus¬ 
sian centered on its true [Fe/H] value with a equal 
to the assigned error. In the recovery, the in¬ 
dividual likelihoods are independent of mass and 
age, as they only depend on [Fe/H]; the likelihoods 
are Gaussians centered on the observed metallic- 
ity with cr given by the individual [Fe/H] mea¬ 
surement errors. Here and in real datasets, the 
impact on computing resources is small, because 
the number of stars in the spectroscopic sample is 
much lower (generally only a few hundreds) than 
the photometric sample, and all the terms needed 


to compute the equivalent of Eq. (13) for the spec¬ 
troscopic PPP are already calculated when solving 
for the photometric PPP. 


To demonstrate the effects of a spectroscopic 
constraint, we show the cases of catalogs Blg and 
Blg:No390; the former is used as template, while 
the latter does not contain measurements for the 
metallicity-sensitive F390W filter. We put our¬ 
selves in the situation where the catalog is gen¬ 
erated with a Gaussian DM distribution, with 0 
mean and a = 0.25 mag, and we use the correct 
prior to marginalize over distance. The results are 
shown in Fig. El In the first and third panels 
from the bottom, we show the results obtained 
without imposing a spectroscopic constraint, while 
in the second and fourth panels from the top the 
constraint is used. The results that include spec¬ 
troscopy are more accurate and precise than the 
results that do not, though they are still not as 
good as the corresponding cases without a dis¬ 
tance spread; this outcome is not surprising, given 
the examples of catalog Burst3 that explored in 
Sect. |7.4[ The improvement is more noticeable 
for catalog Blg:No390 (third and fourth panels 
from the top) than for catalog Blg (first and sec¬ 
ond panels from the top), since the latter includes 
the metallicity-sensitive F390W band. 


8. Summary 

We have introduced a new approach to the 
study of resolved stellar populations via multi¬ 
band photometric observations. The outlined 
framework is based on PPP theory. We solve 
the problem using standard Markov Chain Monte 
Carlo techniques, combined with techniques de¬ 
veloped for medical imaging reconstruction, such 
as sparsity regularization for Poisson Data. 

The underlying idea driving this work was the 
need to simultaneously solve for the IMF slope, 
SFH, and MDF for nearby environments such 
as the Galactic bulge and the Milky Way satel¬ 
lites. We have developed a framework that al¬ 
lows easy inclusion of nuisance parameters, such 
as stellar distance, and demonstrated the impor¬ 
tance of specifying informative priors for these nui¬ 
sance parameters. We have shown how to robustly 
incorporate measurement errors, incompleteness, 
and selection functions within the PPP frame¬ 
work. Our approach is particularly useful when 
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Fig. 11. — Results for catalogs Blg (top two panels) and Blg:No390 (bottom two panels) with a Gaussian 
DM distribution. The first and third panels from the top show the results for the basic recovery. The second 
and fourth panels from the top show the results when a spectroscopic constraint on the MDF is applied. See 
Fig. I^for a description of the individual panels. 
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multi-band data are available; in this case meth¬ 
ods based on CMD gridding can be less straight¬ 
forward to apply. Another advantage of our ap¬ 
proach is the ease with which we can incorporate 
certain types of additional observations, such as 
those coming from independent spectroscopic ob¬ 
servations. 

We have validated our method by simulating 
catalogs with different underlying stellar birth 
functions (number of stars formed per unit mass, 
age, and metallicity) and showing the how well we 
can recover the input values. We have tested the 
outcomes under different assumptions on the pho¬ 
tometric errors, catalog size, selection function, 
available photometric bands, and accuracy of prior 
assumptions on the nuisance parameters. These 
tests demonstrate that our technique recovers the 
input parameters without significant biases, lim¬ 
ited only by the uncertainties in the data. 
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A. The Expectation-Maximization Algorithm 


The Expectation-Maximization Algorithm (EMA) is a general iterative algorithm for finding the 
maximum-likelihood, or in our case maximum-a-posteriori (MAP), parameters of a probability distribu¬ 
tion function when a direct solution is non-trivial. There are many applications of the EMA; we refer the 
reader to Dempster et al. (19771, McLachlan & Krishnan (20081 and Streit (2010) for further details. We 
summarize, without derivation, the EMA steps for PPPs on a discrete space (our grid of models). 

Consider a set of tj, j = 1...M measurements (stellar magnitudes) and a piece-wise constant PPP, i.e. 

KIr{s), where K is the number of grid cells, and Ir is 1 across the r-th cell and 0 elsewhere. 
The likelihood of the data, given A, is: 


M 


K 




1=1 \''=1 


(Al) 


where 



p(t I s) Ir{s) ds 


and 


iTlrl 


Ir{s)ds. 


The idea of the EMA is to introduce a number of latent variables, usually referred to as missing data. In 
this case the missing data are the true values of the stellar parameters, which we indicate with u. The joint 
probability of data and latent variables is given by: 


A) =e 





(A2) 


Combining Eqs. (|Al[) and (A2) it is possible to derive the expression for the probability of the latent 


variables conditional on the data and A: 


p{{u}\{t},X) = 


p({u},{t}| A) 
P(W|A) 


M 

n 


Auj fuj {tj ) 
=^1 Y.r=l Ar fr{tj) 


(A3) 


The E-step consists of taking the logarithm of the joint probability, Eq. (A2), and calculating its expecta¬ 
tion value over the conditional, Eq. (A3). The M-step consists of maximizing the expression obtained in the 


E-step. Without showing it, we note that the M-step requirements lead to an iterative scheme for computing 
Xr until a given convergence tolerance is reached, e.g. until the difference between the log-likelihood at step 


n and n -|- 1 is below a fixed threshold. This algorithm falls in the class of Shepp-Vardi algorithms (Shepp 
&: Vardi] 1982). Further modifications to the EMA are necessary when including the sparsity regulariza¬ 
tion penalty logarithmic function. We will not illustrate those modifications here and refer the reader to 


Lingenfelter et al. (20091, where the recursive updates for the EMA are derived for this particular case. 


B. Catalog Figures 
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Fig. 12. — Same as Fig. but for the BurstS catalog 
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Fig. 13. — Same as Fig. but for the ExP catalog 
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Fig. 14. — Same as Fig. but for the Const catalog 
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